library(tidyverse)
library(sf)
library(leaflet)
library(htmltools)
library(htmlwidgets)
library(raster)
library(gstat)
library(spatial)
library(dplyr)
library(wesanderson)
leaflet_plane <- "+proj=longlat +datum=WGS84"
MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"
traffic <- st_read("traffic.csv", quiet = T) %>%
filter(Lat != "NA")
traffic2 <- st_as_sf(traffic, coords = c("Long", "Lat"), crs = crs(leaflet_plane)) %>%
mutate(total_peak = as.numeric(as.character(Peak_total)))
streets <- st_read("http://bostonopendata-boston.opendata.arcgis.com/datasets/cfd1740c2e4b49389f47a9ce2dd236cc_8.geojson?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D",
quiet = T) %>%
st_transform(crs = leaflet_plane)
south_boston <- st_read("C:/Users/mwill/OneDrive - Harvard University/2020 - Fall/Spatial Analysis/InteractiveMaps/SouthBoston_tract", quiet = T) %>%
st_transform(crs = leaflet_plane)
leaflet(south_boston) %>%
addProviderTiles(providers$Stamen.Toner) %>%
addPolygons() %>%
addPolylines(data = streets)
south_boston_streets <- streets[south_boston,]
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
ggplot(south_boston_streets) +
geom_sf()

south_boston_traffic <- st_join(south_boston, traffic2, join = st_intersects)
south_boston_traffic2 <- south_boston_traffic %>%
group_by(GEOID, NAMELSAD, Associated.Streets) %>%
summarise(total_peak_tract = sum(total_peak))
traf_mean <- mean(south_boston_traffic2$total_peak_tract, na.rm = T)
south_boston_traffic2$total <- ifelse(is.na(south_boston_traffic2$total_peak_tract) == T,
traf_mean,
south_boston_traffic2$total_peak_tract)
south_boston_traffic2$describe <- paste("<p>",south_boston_traffic2$NAMELSAD,
"</p>Intersection:",
str_to_title(
south_boston_traffic2$Associated.Streets),
"<br>",
"Vehicles during Peak Hours:",
prettyNum(south_boston_traffic2$total_peak_tract,
big.mark = ",")) %>%
lapply(htmltools::HTML)
bins <- seq(min(south_boston_traffic2$total_peak_tract, na.rm = T),
max(south_boston_traffic2$total_peak_tract, na.rm = T), by = 1)
pal <- colorNumeric(palette = c("#fee6ce", "#fdae6b", "#e6550d"),
domain = south_boston_traffic2$total_peak_tract)
leaflet(south_boston_traffic2) %>%
addProviderTiles(providers$Stamen.TonerLite) %>%
addPolygons(fillColor = pal(south_boston_traffic2$total_peak_tract),
fillOpacity = 0.4,
stroke = F,
highlightOptions = highlightOptions(fillOpacity = 1),
popup = ~describe) %>%
addLegend(pal = pal,
values = ~total,
title = "Total Vehicles",
bins = 5,
position = "topright")
pal <- colorNumeric(palette = c("#f2f0f7",
"#cbc9e2",
"#9e9ac8",
"#756bb1",
"#54278f"),
domain = traffic2$total_peak)
bins <- seq(min(traffic2$total_peak, na.rm = T),
max(traffic2$total_peak, na.rm = T), by = 1000)
traffic2$describe <- paste("Intersection:",
str_to_title(traffic2$Associated.Streets),
"<br>",
"Vehicles during Peak Hours:",
prettyNum(traffic2$total_peak,
big.mark = ",")) %>%
lapply(htmltools::HTML)
leaflet(traffic2) %>%
addProviderTiles(providers$Stamen.Toner) %>%
addCircles(fillColor = pal(traffic2$total_peak), fillOpacity = .8,
stroke = F, radius = 50,
popup = ~describe,
highlightOptions = highlightOptions(fillOpacity = 1)) %>%
addLegend(pal = pal,
values = ~total_peak,
bins = 5,
position = "topright",
title = "Vehicles Counted during<br>Peak Hours 11/21/19")
traffic_pts_sp <- traffic2 %>%
st_transform(MA_state_plane) %>%
as_Spatial()
tract_poly_sp <- south_boston %>%
st_transform(MA_state_plane) %>%
as_Spatial()
streets_poly_sp <- south_boston_streets %>%
st_transform(MA_state_plane) %>%
as_Spatial()
southie_raster <- raster(tract_poly_sp, res = 10)
gs <- gstat(formula = Peak_total~1, locations=traffic_pts_sp)
idw_interp <- interpolate(southie_raster, gs)
## [inverse distance weighted interpolation]
idw_interp_clip <- mask(idw_interp, streets_poly_sp)
bins <- seq(min(traffic2$total_peak),
max(traffic2$total_peak), by = 1)
pal <- colorNumeric("viridis",
domain = traffic2$total_peak,
na.color = "#00000000")
leaflet(traffic2) %>%
addProviderTiles(providers$Stamen.Toner) %>%
addRasterImage(idw_interp_clip, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal,
values = ~total_peak,
bins = 5,
position = "topright",
title = "Total Vehicles Counted<br>Peak Hours on 11/21/19")